Electron-hole spectra created by adsorption on metals from density-functional theory 
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Non-adiabaticity in adsorption on metal surfaces gives rise to a number of measurable effects, 
such as chemicurrents and exo-electron emission. Here we present a quantitative theory of chemi- 
currents on the basis of ground-state density-functional theory (DFT) calculations of the effective 
electronic potential and the Kohn-Sham band structure. Excitation probabilities are calculated 
both for electron-hole pairs and for electrons and holes separately from first-order time-dependent 
perturbation theory. This is accomplished by evaluating the matrix elements (between Kohn-Sham 
states) of the rate of change of the effective electronic potential between subsequent (static) DFT 
calculations. Our approach is related to the theory of electronic friction, but allows for direct access 
to the excitation spectra. The method is applied to adsorption of atomic hydrogen isotopes on the 
Al(lll) surface. The results are compatible with the available experimental data (for noble metal 
surfaces); in particular, the observed isotope effect in H versus D adsorption is described by the 
present theory. Moreover, the results are in qualitative agreement with computationally elaborate 
calculations of the full dynamics within time-dependent density-functional theory, with the notable 
exception of effects due to the spin dynamics. Being a perturbational approach, the method pro- 
posed here is simple enough to be applied to a wide class of adsorbates and surfaces, while at the 
same time allowing us to extract system-specific information. 



I. INTRODUCTION 



When an atom or molecule interacts with a metal sur- 
face, it generates electronic excitations in the metal. This 
general statement is true because the metal offers a con- 
tinuum of electronic excitations, reaching down to ar- 
bitrary small excitation energies, and hence the adia- 
batic limit for scattering or adsorption of particles at a 
metal surface is ill-definedii Dissipation of some part of 
the particle's energy to the metal surface via electronic 
excitation has been known for a long time, and is the 
physical cause of such phenomena as vibrational damp- 
ing of adsorbate vibrations on metals^ exo-electron emis- 
sion in chemical reactions on metal surfaces,—^ and the 
so-called chemicurrents (see below). Moreover, for the 
sticking of atoms or molecules to a surface, it is nec- 
essary that the adsorption energy is dissipated to the 
substrate. Thus the question arises what role is played 
by electronic excitation processes in the energy dissipa- 
tion during sticking. The ubiquity of (albeit small) non- 
adiabatic contributions to chemical reactions at metal 
surfaces underlines the practical importance of under- 
standing the process of electronic dissipation^ So far, 
most calculations addressing gas-surface interactions as- 
sume that an adiabatic regime would exist, by treating 
the interaction dynamics on a Born-Oppcnhcimer poten- 
tial energy surface. While this assumption may often 
yield plausible results, it is strictly speaking not justified 
for metals, as some electronic excitations must always 
occur, as discussed above. On the basis of calculations 
for dissociative adsorption of hydrogen molecules, it has 
been argued recently that the energy dissipated to the 
electronic system is rather small,— even on a metal, and 
thus has little effect e.g. on the energy dependence of the 
sticking coefficient^ These findings seem to suggest that 



some (probably the larger) part of the adsorption energy 
of the sticking particles is dissipated to excitations of 
the lattice. However, further experimental clarification 
of this point is desirable. The measurement of chemi- 
currents offers an alternative possibility to quantify the 
amount of electronic excitation that takes place during 
chemisorption. In these experiments, the gas-phase par- 
ticles are adsorbed on a metal film, fabricated either by 
previous deposition on a semiconductor substrate or a 
thin insulating barrier layer, forming a Schottky diode^ 
or a metal-insulator-metal (MIM) device^ respectively. 
While particles are being adsorbed, an electric current is 
detected, originating from carriers that have been excited 
sufficiently strongly to overcome the in-built potential 
barrier of the detection device. Typical barrier heights 
are of the order 0.4 - 1 eV. The physics of the underly- 
ing processes is similar to the emission of exo-electrons, 
only differing in the fact that the lower internal barrier, 
e.g. in a Schottky diode, allows for detection of a larger 
fraction of the excited carriers as compared to emission 
into the vacuum. The latter process is possible solely 
for electrons whose energy exceed the work function of 
the metal, typically several eV. The size of the measured 
chemicurrents is dependent on the chemical nature of the 
adsorbing species, and to a lesser extent also on the metal 
surfaced The strong observed isotope effect (so far only 
data for hydrogen versus deuterium adsorption are avail- 
able, differing by about a factor three)^ clearly shows 
that the origin of the measured chemicurrents is related 
to dynamical electronic excitations. The size of the iso- 
tope effect is, so far, the most sensitive test for a theoret- 
ical description of chemicurrents, as it reflects the small 
high-energy tail of the spectrum of electronic excitations. 

In this paper, we arc going to present a theory of chemi- 
currents that allows us to make material-specific prcdic- 
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tions about the spectrum of the electronic excitations 
and the isotope effect. To begin with, one has to keep 
the following observations in mind: While electronic ex- 
citations obviously contribute to the overall energy dissi- 
pation, previous calculation s 6 ' 11 ' 12 have shown that just 
a small fraction of the chemisorption energy is dissipated 
while the particle approaches the chemisorption potential 
well and makes the first round trip. Hence the influence 
of the electronic energy dissipation on the particle's tra- 
jectory is a small (second-order) effect. Moreover, the 
vast majority of the electrons and holes produced by the 
excitation has an excitation energy too small to over- 
come the barrier of the Schottky diode and get detected. 
Hence, the efficiency of an incident particle for creating 
detectable excited carriers in these experiments is tiny, 
less than one percent^^ Therefore, we suggest the ap- 
plication of time-dependent perturbation theory to cal- 
culate these small excitation probabilities, using as input 
the realistic electronic structure of the metal substrate 
obtained from density-functional theory calculations. 

Previously, there have been several other approaches 
in use to calculate chemicurrents: In recent work by 
Lindcnblatt and Pehlke, the dynamics of the full elec- 
tronic system of a slab model of the surface has been 
calculated from time-dependent density functional the- 
ory (TDDFT).— This approach is conceptually very ap- 
pealing, but computationally cumbersome. Moreover, it 
requires high numerical accuracy, as the excitation en- 
ergy of the electronic system, a small number, needs to 
be calculated by comparing huge numbers, the total in- 
ternal energies in the ground state and the excited state. 
In a different approach, the energy dissipated to the elec- 
tronic system has been calculated from the theory of elec- 
tronic friction , 11 ' 13 that was originally developed in the 
context of adsorbatc vibrational damping. However, this 
theory does not yield directly the spectra of the excited 
carriers. Instead, it is assumed that the excitation is 
made up from a superposition of individual electron-hole 
pairs of small excitation energies. These energies are as- 
sumed to be limited to a narrow energy region around 
the Fermi energy of the metal, allowing for a calculation 
of the excitation probability solely based on the prop- 
erties of electronic states right at the Fermi level. For 
higher excitation energies, the probability is calculated 
by referring to the additional assumption of independent 
multiple electron-hole pair excitations that can be de- 
scribed by a Poissonian process^ Thus, the assumptions 
made do not allow one to take into account the electronic 
single-particle properties of the metal under study fur- 
ther away from the Fermi energy. Finally, yet another 
theory has been worked out that allows the calculation 
of chemicurrents by solving the dynamics of an appro- 
priate model system, described by a Newns- Anderson 
Hamiltoniani 15 i 16 ' 17 i 18 ' 19 This scheme is computationally 
lighter than the full TDDFT dynamics, and allows the 
discussion of interesting trends over a variety of model 
systems! 16 ' 17 Material-specific physics can be described 
by matching model parameters to values calculated from 



DFTji£ however, no detailed account of the metallic band 
structure for the excitation of chemicurrents is attempted 
in this theory. 

The time-dependent perturbation theory worked out in 
this paper puts us in position to calculate electron and 
hole excitation spectra separately and directly from the 
Kohn-Sham band structure and wavefunctions obtained 
in DFT. From the spectra, the excitation efficiency per 
impinging particle and its dependence on isotopic mass 
can be calculated without further drastic assumptions. 
The method allows us to make material-specific predic- 
tions while at the same time being computationally af- 
fordable, as only data from ground-state DFT calcula- 
tions are required. The remainder of this paper is orga- 
nized as follows: First, an outline of the theory, including 
a discussion of its limitations, is given. Secondly, results 
are presented for a specific system, H/A1(111), where we 
can compare to the results of full TDDFT calculations 
from the literature. Finally, we discuss the relationship of 
our present theory with other approaches and conclude. 

II. THEORY 

In this section we present the perturbativc approach 
used for calculating the spectrum of electrons and holes 
excited by an adsorbing particle. The motion of the parti- 
cle is described by its classical trajectory. For simplicity, 
our notation assumes a one-dimensional motion along a 
reaction coordinate Q(t), but generalisation to multiple 
dimensions is straightforward. Typically the particle's 
trajectory will make many oscillations in the adsorption 
well before it has dissipated the adsorption energy to the 
substrate. In the following, we consider just the first 
round trip of the trajectory in the adsorption potential 
in the time interval between — r and r. For conceptual 
clarity, we assume that a trajectory unaffected by energy 
dissipation can be used. This assumption is a good first 
approximation in many cases of practical interest, but 
it can be lifted without much difficulty if required. For 
a full description of adsorption, one should decompose 
the oscillating trajectory into many round trips, and use 
a somewhat different energy and corresponding starting 
point for each segment of the trajectory. In our case the 
energy loss during the first round trip is already sufficient 
for the particle to become trapped. 

According to the Rungc-Gross theorem^ time- 
dependent density functional theory (TDDFT) estab- 
lishes a mapping of the electronic many-particle problem 
onto an effective single-particle Hamiltonian of the form 

H°{t) = -^V 2 + V$z{t). (1) 

The effective potential V e e(t) consists of 

V&(t) = V H (t) + V° xc (t) + V ion (t). (2) 

Vi on depends on time implicitly through the particle's 
position Q(t). The charge density needs to be calculated 
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sclf-consistcntly from the solution of the time-dependent 
Kohn-Sham equations. Hence the Hartree and exchange- 
correlation potentials, Vh and V xc , are time-dependent, 
too. 

While TDDFT is in principle exact, we are seeking for 
suitable approximations to arrive at a computationally 
simpler scheme. With this motivation, we assume that 
the true effective potential V e s(t) in the time-dependent 
Kohn-Sham equations can be approximated by the ef- 
fective potential V e s(Q) of a time-independent problem, 
namely, the electronic ground state for a system with the 
particle located at Q(t). In view of the density of ex- 
cited electrons and holes being small compared to the 
overall charge density, this approximation appears plau- 
sible. However, there are cases where deviations from 
the instantaneous ground state are strong; in particu- 
lar those system where exo-emission has been observed. 
From the theory point of view, a class of system can be 
envisaged where the adsorbate, in some points of its tra- 
jectory, must be described (at least in an approximative 
sense) by quantum numbers (such as charge state, spin, 
etc.) that differ from those of the instantaneous ground 
state. For these systems, we anticipate that one may 
need to go beyond the perturbative approach outlined 
here. 

The change of the effective potential due to the scat- 
tering of an incident particle constitutes a localized time- 
dependent perturbation. The chemicurrent consists of 
propagating charge carriers that are detected (in prin- 
ciple) far away from the particle's impact point. These 
carriers are described by the eigenstates \Ei) of the Hamil- 
tonian Hq := H(Qq). In order to describe vibrational 
damping, Qq should be taken as the adsorption height. 
In the case of interest here, we take Qq — * — oo, i.e., Hq 
describes the (semi-infinite) metal with the particle far 
away from the surface. This defines a decomposition of 
the Hamiltonian into an unperturbed part Hq and a per- 
turbation V(Q) = V e e(Q) — Veg(— oo) that depends on 
time implicitly through Q(t). While V(Q) is generally 
not small, its perturbative treatment relies on the fact 
that V(Q(t)) varies rapidly and is different from zero 
only within a short time interval. Here, we are interested 
in electron-hole pair excitations above some threshold en- 
ergy hvs > 0.5 eV due to the high-frequency components 
of the perturbation. For the validity of the perturbative 
treatment, it is required that the response of the elec- 
tronic system on the time scale set by Vg 1 is sufficiently 
weak. Applying first-order time-dependent perturbation 
theory, the transition amplitude for an excitation of an 
electron from an occupied state i into an unoccupied state 
j is given by 



Pij (t) = (ej \V(Q(t))\ Ei ) exp(i(Ej - Si)t/h), 
and the excitation spectrum is obtained from 
f T dt 
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Note that uj is an energy throughout this paper. For a full 
round trip of the (dissipationless) trajectory, V vanishes 
both in the initial and final (again unperturbed) state, 
and we may integrate by parts to obtain an expression 
for the spectrum of electron-hole pairs, 
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exp(i(sj — ei)t/K). (6) 



Extending our theory to finite electronic temperature by 
introducing Fermi occupation factors /(e), one has to 
consider the possibility of transitions where the metal 
electrons loose energy by going from occupied states 
above £p to unoccupied states below £p. Thus, we need 
to define P ex (w) = for uj < 0. Following ref. 0, eq. ([5]) 
needs to be generalized to 
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X5(u>- (Sj -£i))6{u). 

Interpreting defect electrons below the Fermi energy £p 
as holes, first-order perturbation theory gives explicit ex- 
pressions for the electron and hole spectra 
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(8) 



The total energy loss of the trajectory is defined by 

AE = dujLuP cx (u). (9) 
Jo 

In the following section, we will discuss the spectra for 
a half round trip in the adsorption well rather than a 
full one, in order to compare to previous results in the 
literature. For a dissipationless trajectory, the ways 'in' 
and 'out' are equivalent by time inversion symmetry, and 
within our approach we can take the spectrum of the half 
round trip just being equal to one half of the spectrum 
for the full round trip. Also, we note that time-reversal 
symmetry implies Ay = as |e, — Ej\ — > 0. But even 
without this assumption, the expression eq. ([6]) remains 
mathematically well-defined for a half round trip, as the 
matrix element in the integrand vanishes at the turning 
point of the trajectory. For evaluating Ay, we evaluate 
the integrand at discrete timesteps, and use 
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AQ 



dQ 

-exp(z( £ , 



Ei)t/h). 

(10) 
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The DFT calculations are carried out at finite electronic 
temperature T, and the same temperature is used in the 
Fermi functions in eq. and (JSj) . In order to plot 
the calculated spectra, it is necessary to replace the 5- 
functions in eq.s ((7|) and ([8]) by a function of finite width. 
In practice, we replace the 5-function by —df/de. The 
same electronic temperature T as in the DFT calcula- 
tions is used for the broadening. Matrix elements with 
\si — £j\ < k-gT have to be handled with care. Although 
the limit \£{ — £j\ — > is mathematically well-defined, the 
numerical results are not reliable due to rounding errors. 
In the plots, and also in the calculations, this energy re- 
gion is therefore omitted. The total energy loss is not 
affected by this omission, since the contribution of such 
small-energy transitions to AE is negligible. (For a full 
round trip, P cx (0) = 0). 

We note that the approach presented here can be ex- 
tended beyond first-order perturbation theory. This is 
possible if the electron-hole pairs excited by the adsorb- 
ing particle can be treated as independent bosons. In 
this case, the electron-hole spectrum (but not the sep- 
arate spectra for electrons and holes) can be calculated 
analytically. The method is applicable to weak and swift 
perturbations^ An extension to strong perturbations 
slowly varying in time has been suggested as well, where 
the wavefunctions used to calculate the matrix clement 
in eq. (|10p are replaced by the ground-state wavefunc- 
tions calculated at each Q-valuei 22 ' 23 For the system un- 
der study, we find that first-order perturbation theory 
is sufficient (see Section IV). The theory presented here 
has some similarities with, but also important differences 
to the theory of electronic friction worked out earlier^ 
To see this, one assumes the separability of the matrix 
elements appearing in eq. © (see ref. [Tj for a discus- 
sion of this point), and takes the quasi-static limit. This 
limit implies that the excitation frequencies of interest 
are sufficiently small, such that the states |e-^ ) and (£j\ 
are essentially indistinguishable from states at the Fermi 
surface, \ef)- The resulting theory is formally identical 
to the theory of electronic friction, with the important 
difference that this theory uses wavefunctions changing 
adiabatically with the position of the adsorbate, while 
we use wavefunctions of the bare substrate to calculate 
the matrix elements. We remind the reader that the the- 
ory of electronic friction starts from an expression for the 
imaginary part of dynamic self-energy A(Q,lo), 

ImA(Q, W ) = (e s — Si j 

x (/(*) -ffa)) Siw-te-etMll) 

see ref. 0. Note that this definition is somewhat analo- 
gous (although not equivalent) to the definitions in cq. ([6|) 
and ||7J). Whereas the present definition of P ex uses the 
squared modulus of a time integral, the definition of 
ImA starts from the square moduli of matrix elements. 
Hence a friction coefficient rj(Q) can be defined for any 
point Q(t) along the trajectory by considering the limit 



lim^^o mA(Q, uj)/u and replacing the states |ej) and (ej\ 
by the states at the Fermi energy. Taking into account 
the proper normalization of these states on the energy 
shell one can show that this limit exists, i.e., the fric- 
tion coefficient is mathematically well-defined. A similar 
reasoning employing the quasi-static limit can be used 
to establish that the total energy loss, eq. ©, is a well- 
defined quantity also in our theory, both for a full and a 
half round trip. The limit of P ex (uj) for u> — > can be 
shown to exist, as the expression can be decomposed into 
the two factors, |Ay \ 2 /\£j-£i\ and (/(e i )-/(e i ))/|ej-e i |, 
which both have a finite limit for \ei — £j\ — * 0. For the 
first factor, the proof is analogous to the derivation of 
the finiteness r\ in electronic friction theory. The second 
factor is finite for any non-zero electronic temperature 
T. When considering the total energy loss, one may even 
send T — > 0; the combined limiting process ui — > 0; T — > 
for the integrand luP cx (lu) remains well-defined for any 
fixed ratio u)/(k^T). 

In extension over electronic friction, the present theory 
allows to us address the excitations of a single electron- 
hole pair at any finite energy, even further away from the 
Fermi level. Compared to adsorbate vibrational damp- 
ing, chemisorption is accompanied by stronger variations 
of the effective electronic potential, and hence compo- 
nents of higher frequency, also much higher than the vi- 
brational frequency in the adsorption well, may become 
excited non-resonantly. The present theory is applicable 
under the following conditions: The energy loss AE dur- 
ing one round trip must be a small fraction of the total 
energy of the adsorbing particle, measured from the bot- 
tom of the adsorption well, to justify the 'unperturbed 
trajectory' approximation. Moreover, the particle is re- 
quired to move sufficiently rapidly to justify the use of 
static (rather than adiabatic) electronic wavefunctions 
in evaluating the matrix elements. This means that the 
excitation frequency of interest, v = \ej — £i\/h, should 
fulfill the condition v ^> (2ttt) , where r is the duration 
of the impact (in scattering) or of half a round trip (in 
adsorption) . The latter condition is fulfilled for electronic 
excitations detected by a Schottky diode (hv > 0.5 eV, 
i.e. hvs = 0.5 eV) even for light particles such as thermal 
hydrogen atoms, where t ~ 20 fs. 

Hence, the present theory seems to be more appropri- 
ate for the high-energy part of the excitation spectrum, 
detectable as chemicurrents, than the electronic friction 
formalism. Moreover, the present, extended theory has 
some technical advantage concerning the so-called spin 
transition: In most practical implementations of DFT, 
the electronic correlation effects between the electrons at 
the adsorbing particle, e.g. hydrogen, and the metal elec- 
trons are described in an approximate way by solutions 
that break spin symmetry, by assigning a spin moment to 
the free hydrogen atom. During adsorption, spin symme- 
try is restored. In such ground-state DFT calculations, 
the spin transition has a formal analogy to a first-order 
phase transition, and thus leads to a ^/Q-type depen- 
dence of the effective potential near the transition point 



on the reaction path. In the theory of electronic friction, 
this singularity results in a divergence of A(Q, oS), render- 
ing the theory inapplicable.— In the present theory, the 
Ay, and hence P cx (lo), remain well-defined, since the in- 
tegrand in eq. (fit))) diverges at most like Q -1 / 2 , and hence 
remains intcgrable. Phrased in other words, the abrupt 
change of the spin polarization in a ground-state DFT 
calculation and the quasi-static limit are incompatible. 
The present theory circumvents this difficulty. However, 
we mention that in a time-dependent treatment of spin, 
both in TDDFT- 2 or in the mean-field approximation 
to the Newns- Anderson model,— the spin polarization 
changes smoothly, and this difficulty doesn't arise any- 
way. 



III. CALCULATIONS 

The above concepts are illustrated by calculations for 
an H atom impinging on the on-top site on an Al(lll) 
surface. We chose this system because we can compare to 
recent results where the full dynamics has been treated 
within TDDFT^ Here, we perform static DFT calcula- 
tions within the generalized gradient spin density approx- 
imation (GGSDA). We employ the PBE-GGA exchange- 
correlation functional 2 ^ together with norm-conserving 
pscudopotcntials. We model the system by an Al(lll) 
slab of 12 layers with a 2v / 3 x 2-\/3 unit cell. The slab is 
relaxed with respect to the interlayer distances, while we 
keep the lowest 3 layers unrelaxed. Electronic wavefunc- 
tions are calculated using the pwscf-code~ for 35 po- 
sitions of an H atom above the surface on the on-top 
position of one of the substrate aluminium atoms, and 
for the substrate slab alone. The cut-off energy for the 
plane-wave expansion of the wavefunctions is chosen to 
be 21 Ry, and a 4 x 4 x 1 Monkhorst-Pack k-point mesh 
is used for the Brillouin zone integration. The cut-off en- 
ergy was tested by calculating the chemisorption energy 
of H for the on-top position on Al(f 11), once using the 
cutoff of 21 Ry, yielding 1.83 eV, and once using a cutoff 
of 40 Ry, yielding a chemisorption energy of 1.86 eV. The 
difference of 30 meV is acceptable in our opinion, when 
compared to the much larger calculation times for the 
higher cutoff. TDDFT calculations did not even reach 
this cutoff if a reasonable calculation time is a necessity, 
but the results from static calculations compare well to 
our numbers (1.90 vs. 1.87 eV for 40 and 20 Ry cutoff, 
respectively)^ A smearing of the occupation factors of 
1 mRy using a simple Fermi-Dirac broadening scheme 
is used, which corresponds roughly to the experimental 
temperatures^ Using the aformentioned parameters we 
obtained a lattice constant of 4.061 A for bulk Al, while 
the experimental value is 4.050 A. 

Fig. [T] shows the Kohn-Sham band structure of an 
Al(lll) slab with a (1 x 1) unit cell resulting from our cal- 
culations. Comparison with the literature^ shows good 
agreement. The lowest state at the T-point is about 11 eV 
below, the band bottom at the ill-point about 5 eV, and 
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FIG. 1: Kohn-Sham band structure of an Al(lll) slab con- 
sisting of 12 layers in a (1 x 1) unit cell calculated within 
PBE-GGA. A Monkhorst-Pack k-point mesh of 16 x 16 x 1 
points has been used to calculate the charge density and the 
self-consistent potential. The Fermi energy is set to 0. 



at the iGpoint about 3 eV below the Fermi energy. 

Fig. [5] shows the chemisorption potential of hydrogen 
above the on-top site of Al(lll). The ground-state po- 
tential energy surface (PES) consists of two branches that 
intersect at a distance of about 2.6 A above the sur- 
face. For larger distances, the DFT ground state is spin- 
polarized, while the electronic system is non-polarized 
when the H atom comes closer. As seen in fig. O the 
electronic spin polarization at the H atom vanishes be- 
low this point. 

Fig. [3] shows the result of a population analysis for 
the atomic Is orbital of the H atom using the Lowdin 
approach implemented in pwscf. The DFT ground-state 
wavefunctions of either spin for each position along the 
trajectory were projected onto the H Is orbital. Below 
the spin transition point at about 2.6 A both the spin- 
up and the the spin-down orbital shows a population of 
about 0.6 electrons. Further away from the surface, the 
H atom is spin-polarized, i.e., the spin- up orbital gets 
fully populated while the spin-down orbital gets depleted. 
Right above the transition point the spin polarization 
shows a sharp, square-root-like increase. 

The classical trajectory for vertical impact of the 
H-atom onto the on-top site of Al(lll) is calculated 
from Newton's equations for an initial kinetic energy of 
60 meV, and starts 3.16 A above the surface. Thereby, 
the small effect of electronic friction on the trajectory is 
neglected. The total time needed for one half round trip 
is 20 fs. 

The numerical evaluation of eq. p0[) proceeds as fol- 
lows: The derivative dV/dQ is calculated numerically, 
using the first derivative of a third-order spline inter- 
polation. The potential used is the effective potential 
from the ground-state DFT calculations. The matrix el- 
ements in eq. (|10[) arc computed in two steps: First we 
perform a summation in real space to calculate AV^e^), 
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FIG. 2: (Color online) Chemisorption potential of a H atom 
above the on-top position of an Al(lll) surface. The spin 
polarization in the ground-state DFT calculations vanishes 
at about 2.6 A. The dashed lines are two Morse potentials 
fitted to either branch of the PES. 



then we make a Fast Fourier Transform to wave-vector 
space and multiply by (sj \ . The wavefunctions are those 
of the initial state at Q(— oo), but omitting adsorbate 
wavefunctions. Likewise, the energies in the exponen- 
tial function in eq. ([TU)) are the Kohn-Sham energies of 
the initial state. The time integration is then done us- 
ing a spline interpolation of the matrix elements. For 
the technical implementation the state indices i and j 
incorporate band index, k-point and spin. Only transi- 
tions within the same k-point and spin channel are con- 
sidered. However, the restriction regarding the k-points 
can be lifted^ A weight factor for the summation over 
k-points also has to be introduced. 



IV. RESULTS 

The main result of the theory presented here are the 
spectra of the electrons and holes excited by a hydro- 
gen atom impinging with thermal kinetic energy. From 
the part of the spectra exceeding the Schottky barrier, 
we directly derive the predicted yield of the measured 
chemicurrent. In order to study the isotope effect system- 
atically, we not only carry out calculations for H and D, 
but also for fictitious atoms of mass 0.1m# and 0.01m# . 
The spectra are shown for a half round trip in the adsorp- 
tion potential (until the particle has reached the turning 
point on the repulsive branch of the PES), to be able to 
compare with earlier results.— We note that in the ex- 
periments the particles perform many round trips before 
they have dissipated the chemisorption energy and come 
to rest. Hence the average number of electron-hole pairs 
per incident hydrogen atom in an actual experiment will 
be significantly higher than the values displayed in the 
figures. 
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FIG. 3: (Color online) Ground-state occupation of the spin- 
up and spin-down orbital of an approaching H atom vs. its 
distance from an Al(lll) surface. Far away from the surface, 
the electronic charge at the H atom, corresponding to one 
electron, is fully spin-polarized. Both spin-up and spin-down 
orbitals become equally populated after the spin transition 
point. The spin polarisation shows a \/Q-behavior close to 
the transition point. 

Fig. [4] shows the spectra of electrons (positive ener- 
gies) and holes (negative energies) for a half round trip 
of an H atom. We notice that the spectra of electrons 
and holes are quite similar. Moreover, differences in the 
spectra of spin-up electrons (same spin polarization as 
the impinging H atom) and spin-down electrons (oppo- 
site spin polarization as the impinging H atom) are only 
marginally different when plotted on a logarithmic scale. 
The same holds for the holes of either spin. The full and 
dashed lines in fig. [3] indicate the results of convergence 
tests. While a metal has a continuous spectrum of excita- 
tions, the slab used in the DFT calculations has a large 
but finite number of states, due to the finite thickness 
and the number of k-points used to sample the Brillouin 
zone. For convergence tests, we increased this sampling 
from a 4 x 4 x 1 k-point mesh to a 6 x 6 x 1 mesh. (This 
refers to the Kohn-Sham wavefunctions \ei) used to cal- 
culate the matrix elements in eq. (p~0|) ; the potential is 
assumed to be converged already for the smaller k-point 
set.) Fig. 2] shows that the spectra are almost identical, 
albeit the 6x6x1 calculations gives somewhat smoother 
curves. The energy loss AE is 23 meV for the smaller 
k-point set and 20 meV for the larger one. The error ap- 
pears to be tolerable, and for the following calculations, 
the smaller set was used. Note that the energy loss is 
just a small fraction, about 1%, of the chemisorption en- 
ergy, see fig. [2] Hence, the use of an elastic trajectory 
in calculation of the matrix elements appears to be justi- 
fied. Fig.[5]shows the spectrum of the electron-hole pairs. 
Here we compare the results of first-order perturbation 
theory, eq. ([7]), to the treatment by Muller-Hartmann 
and co-workers^ that sums up all orders in the pertur- 
bation expansion. We find that the two approaches yield 
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FIG. 4: (Color online) Spin-resolved excitation spectra of elec- 
trons (positive energies) and holes (negative energies) for an 
impinging H atom with an initial kinetic energy of 60 meV on 
the on-top position of an Al(lll) surface after one half round 
trip in the chemisorption well. Values close to the Fermi en- 
ergy (|e— ef\ < 0.013 eV) are not plotted because of numerical 
inaccuracies in this energy range. 



very similar spectra when plotted on a logarithmic scale. 
This demonstrates that first-order perturbation theory is 
sufficient for the high-lying excitations with small excita- 
tion probability, which are of interest to the calculation 
of chemicurrents. Another remark is of interest here: If 
only the combined electron-hole spectrum were known, 
one would be tempted to assume that electrons and holes 
are distributed symmetrically around the Fermi energy. 
Comparing to fig. [H we find that this assumption is rea- 
sonable in the present case (but not for the stronger ex- 
citations by fictitious atoms lighter than H, see below). 

Next, we investigate the isotope effect. The electron 
and hole spectra for deuterium are shown in fig. [SJ They 
are qualitatively similar to those of H, but the excita- 
tion probability decays more strongly with increasing ex- 
citation energy of the charge carriers. Inspired by the 
work of Lindcnblatt and Pehlkep£ we also used fictitious 
isotopes of smaller mass, 0.1m# and 0.01m#, in order 
to get stronger excitations that facilitate the analysis. 
The spectra are collected in fig. [7J Similar to previ- 
ous work;2& we fit an 'effective temperature' T e g to the 
high-energy tails of the spectra, using a two-parameter 
fit P(e) = Aexp(— e/fcsTeff)- For the hole spectra, the 
fitted regions are shown in fig. [7J The temperatures ob- 
tained from the fit are given in table U 

The quite large error bars were estimated by varying 
the fitting window. They indicate that the distributions 
are not truly exponential. A similar uncertainty of the fit- 
ted temperatures was also noticed in previous worki 18 i 28 
We note that the temperatures fitted for the excited elec- 
trons are systematically higher than for the holes. This is 
particularly obvious for the light isotopes, but also holds 
for H and D, despite the nearly symmetric appearance 
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FIG. 5: (Color online) Spin-resolved excitation spectrum of 
electron-hole pairs for an impinging H atom with an initial 
kinetic energy of 60 meV on the on-top position of an Al(lll) 
surface after one half round trip in the chemisorption well. 
The curves calculated in first-order perturbation theory are 
in good agreement with those obtained by a resummation of 
the perturbation expansion following Miiller-Hartmann et al. 
(MHRT). Values close to zero energy (e < 0.013 eV) are not 
plotted (and omitted in the calculations) because of numerical 
inaccuracies in this energy range. 



Isotope 


Perturbation theory 
T c e ff [K] T': s [K] 


TDDFT22 

T: s [K] T c h ff [K] 


Deuterium 
Hydrogen 

#Exp 

H : m — O.lrriH 
H : m — 0.01m H 


iioo-ll looo-fa 

1400 ~ 300 1300 ~ 300 

lwu +700 ±JUU +700 

3000^ 2600^ 
8900-§2f 7300-Hf 


yuu +ioo 1±UU +100 

2700^200 2 3 00-ia 
9300^ 7800-f£ 



TABLE I: Slope parameter T e g of the electron and hole spec- 
tra for different isotope masses. The middle column contains 
the fits to our data, and the right column contains fits to 
data from TDDFT calculations.— The value -Hex P is the ex- 
perimental value for H adsorption on Ag. The error bars for 
the fits are quite large. The reason is that the fits can be 
made over a wide energy area, starting and ending on differ- 
ent points, and thus are rather arbitrary in a certain range. 



of the spectra in fig. [5] Within the error bars, the tem- 
peratures obtained from perturbation theory are in rea- 
sonable agreement both with experiment and with the 
TDDFT calculations. Although the effective tempera- 
tures are helpful for identifying trends, an analysis of the 
isotope effect on this basis is not conclusive due to the 
large error bars, and it is safer to study the isotope effect 
in the total energy loss, and in the integrated probabil- 
ity of excited carriers above the Schottky barrier. Fig. [5] 
shows that the linear relationship between the total en- 
ergy loss and the mass scaling parameter (m/m^-) -1 / 2 . 
This can be rationalized in view of electronic friction the- 
ory. In this theory, the friction coefficient depends only 
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FIG. 6: (Color online) Spin-resolved excitation spectra of elec- 
trons and holes for D with an initial kinetic energy of 60 meV 
after one half round trip in the chemisorption well. The spec- 
trum is qualitatively similar to the spectrum for H, shown in 
fig- HI but the decay is much faster. 

on the electronic properties of the system, and hence 
would be the same for all isotopes of the same chemical 
element. The frictional force is thus proportional to the 
velocity of the impinging particle. By plotting the energy 
loss as a function of (m/m^) -1 / 2 , one would therefore 
expect all data points to fall on a straight line. This 
also holds in the present, more elaborate theory. The 
integrated probabilities N e and Nh of excited carriers 
above the Schottky barrier are shown in tabic [III where 
a barrier height of 0.48 eV has been assumed, typical of 
experimental Schottky barriers for Ag on Sii 10 ' 29 



Isotope 


Electrons 


Holes 


Deuterium 


0.000246 


0.000175 


Hydrogen 
N/Nn 

N/ Nr>, Experiment 


0.000781 
3.17 


0.000561 
3.21 
3.7 ± 0.7-iSi^ 


H : m — O.lmg 

N/Nh 


0.021502 
27.53 


0.014713 
26.23 


H :m = O.OlrriH 
N/Nh 


0.088775 
113.67 


0.064376 
114.75 



TABLE II: Average number of electrons and holes created per 
impinging H isotope with an initial kinetic energy of 60 meV, 
per half round trip. 

In electronic friction theory, the excitation probabil- 
ities are assumed to be Poissonian distributions^. All 
moments of such a distribution a solely determined by 
its first moment. It follows that all excitation spectra (of 
different isotopes) for a given chemical species impinging 
on a given substrate are described by a one-parameter 
family of curves. Then, a single parameter, e.g. the 
energy loss AE, is sufficient to characterize the spectra 
completely. Here, we test if such a simple description 



FIG. 7: (Color online) Spin-resolved excitation spectra of elec- 
trons and holes for different impinging isotopes of H with an 
initial kinetic energy of 60 meV after one half round trip in the 
chemisorption well. Spectra for three different masses with 
m = ran, rn — O.lmn and m = 0.01m_H are shown, where 
lower masses give spectra with a larger amplitude. Decaying 
exponentials with A exp(—e/kBT e g) are fitted to the hole part 
of the spectra. There is some freedom in choosing the energy 
window for the fit, and hence the slopes are associated with 
large error bars. The lines shown in the plot correspond to 
T eff = 1300 K, 2600 K, and 7300 K, respectively. 



of the spectra also applies to the results of our theory. 
The scaling properties of AE have already been estab- 
lished above. If one supposes a proportionality of the 
effective temperatures T Q g of the high-energy tail and 
the total energy loss AE, one would expect that a loga- 
rithmic plot of N e or Nh shows a linear dependence on 
(m/m//) 1 / 2 . Figure [9] demonstrates that this is a good 
approximation for masses between 2tojj and 0.01m# , i.e., 
it is fair to approximate the energy spectra by a single- 
parameter family of curves. The ratio of the integrated 
probabilities N e / h for H and D atoms is about 3.2. (3.17 
if calculated for the excited electrons, 3.21 for the holes, 
see table [TTJ) . Experimentally, the hole excitation prob- 
abilities are more accessible, since Schottky barriers on 
p-type silicon can be fabricated with a high quality of 
the metal-semiconductor interface. The experimental ra- 
tio of the yields for H and D for Ag/Si Schottky diodes 
is 3.7 ± 0.7 ) 10 ' 26 which agrees well with our results. 



V. DISCUSSION 

At present, the possibilities to test our results against 
experiment are rather limited. Quantitative comparison 
is possible so far only for the isotope effect: Here, our 
results for Al are in good agreement with the measure- 
ments using Ag/Si-Schottky diode o 10 ' 26 and MIM sen- 
sors&. Some experimental information is also available 
about the slopes of the excitation spectra from measure- 
ments with Schottky diodes made from the same mate- 
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FIG. 8: (Color online) Total energy loss vs. the mass scaling 
parameter 1 / ^Jm/rriH- rn denotes the mass of the adsorbing 
atom. The line is a linear fit to the data. 



FIG. 9: (Color online) Yield of the electron (red, upper curve) 
and hole (blue, lower curve) excitation spectra plotted versus 
the mass-scaling parameter \J m/mn- rn denotes the mass of 
the adsorbing atom. The lines are linear fits to the data. 



rial but with slightly different /- ^-characteristics. The 
temperatures estimated from these measurements are in 
reasonable agreement with the slopes of our calculated 
spectra. However, we are quite sceptical whether a sin- 
gle temperature would still be a good description if one 
had better experimental access to the spectrum. 

A more detailed comparison of our results is possible 
with previous TDDFT calculations for the same system, 
H on Al(lll). In non-spinpolarized calculations, a to- 
tal energy loss of 30 meV and 40 meV was found after 
a half and a full round trip, respectively^ This is in 
reasonable agreement with the values of 23 meV and 46 
meV calculated in our approach. In addition, the lin- 
ear dependence of the energy loss on the mass scaling 
parameter (cf. fig. [8]), and the effective temperatures 
of the spectra (cf . table [TJ) obtained in both theories are 
consistent. However, more recent spin-polarized TDDFT 
calculations yielded energy losses of about 100 meV and 
185 meV for the half and full round trip, respectively.— 
These findings from the TDDFT calculations suggest 
that the spin transition is quite important for the en- 
ergy loss. While our theory in principle includes the spin 
transition, it seems to underestimate its role in the en- 
ergy dissipation. Preliminary calculations where we use 
a smoothened and delayed spin transition (similar to the 
findings of ref . Il2h instead of an abrupt one show a higher 
energy loss in perturbation theory, too. Moreover, in 
the TDDFT calculations, it was assumed for technical 
reasons that the total spin polarization (integrated over 
both the hydrogen atom and the substrate) is conserved, 
while in our ground-state DFT calculations the total spin 
polarization changes from about one electron for the hy- 
drogen far away from the surface to zero after hydrogen 
adsorption. This methodological difference might be the 
reason why Lindcnblatt and Pchlkc find a shift (in en- 
ergy) between the spin-up and spin-down spectral after 
a half round trip, while our spectra are almost identi- 



cal for both spins. However, this discrepancy does not 
arise after a full round trip, where spin-dependent fea- 
tures are absent in both methods. We may also com- 
pare our results to earlier work on H/Cu(lll) using the 
electronic friction formalism. 1:L i 13 i 14 This theory assumes 
that the energy loss proceeds by multiple excitation of 
very small energy quanta, with the excitation energy of 
a single quantum ultimately tending to zero, while the 
total energy loss remains finite in this limit. Under this 
assumption, the energy loss due to friction can be de- 
scribed by an instantaneous friction force in each moment 
of time, and the materials properties of the substrate only 
need to be known in vicinity of the Fermi energy in or- 
der to calculate the friction coefficient. This is different 
from our approach, where we calculate the probability to 
deposit a finite amount of energy to the substrate within 
a single excitation. Both approaches are complementary: 
While our method focusses on the spectral properties of 
the excitations, including electronic band structure de- 
tails above or below the Fermi energy, our formulae are 
applicable only to a half or full round trip of the parti- 
cle, i.e. at turning points of the trajectory. The electronic 
friction formalism, on the other hand, allows one to break 
down the energy loss (at least within the limitations im- 
plied in this formalism) to every point along the parti- 
cle's trajectory. In order to obtain the excitation spectra 
from the friction approach, additional assumptions need 
to be made, such as the separability of matrix elements, 
and the statistical independence of multiple excitations. 
The spectra (for H/Cu(lll)) obtained after making these 
approximations are essentially exponential distributions. 
Our results from first-order perturbation theory are qual- 
itatively similar, but cannot be fit by a single exponen- 
tial. This implies that the picture of multiple excitations 
made up from energy quanta close to zero may be an 
oversimplification. We note that a connection between 
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both approaches can be made on the basis of the over- 
all energy loss per round trip they predict. From the 
more detailed spectra obtained from our approach, it is 
in principle possible to determine a posteriori an effective 
friction coefficient that would lead to the same amount of 
energy loss as obtained by integrating over our excitation 
spectra. This effective friction coefficient could then be 
used to re-run the trajectory to obtain new spectra, and 
the procedure could be iterated until consistency between 
the energy loss of the trajectory and the loss implied by 
the spectra has been reached. 

Recently, non-adiabatic effects in adsorption 
have been studied in the framework of a time- 
dependent Newns- Anderson model treated in mean-field 
approximation . 15 i 16 ' 17 i 18 ' 19 While this model is originally 
motivated by the surface chemistry of transition metals 
and noble metals, it may still make sense to compare 
the results to our calculations for H on Al(lll). The 
model requires a fit to the data supplied by DFT 
calculations, i.e. t he p rojected density of states onto the 
H orbitals. In ref. [l8|, Miziclinski et al. obtained spectra 
for H/A1(111) that show clear differences between 
electrons and holes, and between both spin orientations. 
However, the results are rather sensitive to the model 
parameters. One can observe that the parametrization 
of the values needed for the model calculations does not 
fit perfectly, especially far away from the surface where 
ionization and affinity level of H are about equally far 
below or above the Fermi energy of Al^ For this case, 
Miziclinski et al. predict electron and hole spectra that 
are approximately equal, and which have no important 
spin characteristics, cf. ref. [lj^ for metals other than 
Aluminium. We think it can be justified to argue that 
H on Al(lll) also corresponds in the language of this 
model to such a system, so that one would expect 
electron and hole spectra that are approximately equal. 
This is found in our calculations, too. However, their 
model also predicts that electron and hole spectra may 
be strongly different for other combinations of substrate 
and adsorbate. In order to probe the different regions 
of the parameter space of the Newns-Anderson model 
discussed in ref. Ilq . i.e., strongly electronegative or 
electropositive adsorbates, additional DFT studies for 



such material systems would be desirable. 



VI. CONCLUSIONS 

In this work, we implemented first-order time- 
dependent perturbation theory to calculate excitation 
spectra of electrons and holes in a metal substrate due 
to non-adiabatic effects in adsorption, starting from the 
ground-state electronic potential and wavefunctions ob- 
tained from density functional theory. The calculated 
results for adsorption of hydrogen on aluminum are con- 
sistent with the experimental data obtained for hydrogen 
on noble metals. In particular, the present theory repro- 
duces the observed isotope effect. In contrast to the elec- 
tronic friction formalism, our approach includes excita- 
tions over finite energy differences, and thus takes the ef- 
fects of single-particle band structure into account in the 
excitation spectra. Our results are in agreement with re- 
sults from the (computationally much more demanding) 
TDDFT method, with the notable exception of spin ef- 
fects. On the other hand, our method using perturbation 
theory is much faster and affordable for a wider range of 
systems, as it requires only ground-state DFT calcula- 
tions. These calculations are widely used in the scientific 
community, e.g. in the calculation of potential energy 
surfaces within the Born-Oppcnhcimcr approximation. 
The method presented here can be implemented in exist- 
ing (static) DFT codes as a post-processing tool, and is 
applicable in cases where the full dynamics of the many- 
particle system remains close to the Born-Oppcnhcimcr 
surface. 
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